home *** CD-ROM | disk | FTP | other *** search
/ BMUG PD-ROM 1996 Fall / BMUG Fall'96 PD-ROM.iso / Education / Math / MathPad 2.4 / Examples / minimize < prev    next >
Text File  |  1996-04-10  |  3KB  |  104 lines
  1. -- Perform multidimensional function minimization.
  2. ~ The user must supply:
  3.  
  4. "f(parms)"  The function to be minimized where
  5.             "parms" is an array of parameters.
  6. "guess"     An array  of initial parameter values.
  7. "deltas"    Initial step size for the parameters.
  8.             If the parameters have different
  9.             scale sizes, "deltas" should define an
  10.             array with a step size for each,
  11.             otherwise a scalar can be used.
  12. "tol"       The termination criteria based on the
  13.             amount of change in parameter values.
  14.             A tol of 1.0e-6 gives about 3 digit
  15.             accuracy.
  16.  
  17. The routine "minimize" steps until a minimum is found.
  18. After minimize is run the following globals are set:
  19.  
  20. "minp"  The solution parameter array
  21. "p"     The final simplex vertices
  22. "y"     The function values at each vertex
  23.  
  24. ~
  25. ------------ Example --------------------
  26. -- fit a curve to data by minimizing the sum of the square of the differences between data and fit.
  27.  
  28. fit(a,x) = a[1]+a[2]*cos(a[3]*x)
  29.  
  30. data=read(xydata); xx[i]=data[i,1]; yy[i]=data[i,2] 
  31.  
  32. err(a) = sum((fit(a,xx[j])-yy[j])^2,j,1,count(data))
  33.  
  34. -- set up globals for minimize
  35. f(a) = err(a)
  36. guess = {80,-75,1}
  37. deltas = .5
  38. tol = 1e-5
  39.  
  40. plot data
  41.  
  42. minimize:;    minp:{80.2,-73.4,0.895}
  43. plot fit(minp,X)
  44.  
  45. ----------- minimization routines -------
  46. -- These routines use the "downhill simplex method". A simplex is a N dimensional geometrical figure with N+1 vertices. The algorithim evaluates the function at each vertex and steps the simplex toward the minimum by: 1) reflecting a vertex away from the function's highest point. 2) reflection and expansion away from the highest point. 3) contraction away from the highest point. 4) contraction in all dimensions toward the lowest point.
  47.  
  48. -- step until parameters change by less than tol
  49. minimize = init,
  50.            step,
  51.            step while pchange > tol,
  52.            minp:=psum/m
  53.  
  54. pchange = sum((abs(p[lowest]-p[highest])/
  55.            (abs(p[lowest])+abs(p[highest])))[ii],ii,
  56.               1,n)
  57.  
  58. -- move simplex one step
  59. step =  rank,
  60.  try(-alpha),                -- reflect from highest
  61.  try(gamma) when ytry ≤ y[lowest], -- keep going
  62.  (ysave:=y[highest],
  63.   try(beta),                 -- contract from highest
  64.   shrink when ytry ≥ ysave) when ytry ≥ y[high]
  65.  
  66. -- find lowest, highest and next highest vertex
  67. rank = lowest:=1,high:=1,highest:=2,i:=1,
  68.    (lowest:=i when y[i]<y[lowest],,
  69.     (high:=highest,highest:=i) when y[i]>y[highest],,
  70.     high:=i when y[i]>y[high] and i ≠ highest,,
  71.     i:=i+1) while i≤m
  72.  
  73. -- move highest vertex through face by fac
  74. try(fac) = ptry:=psum*(1-fac)/n-
  75.                   p[highest]*((1-fac)/n-fac),
  76.     ytry:=f(ptry),
  77.     (y[highest]:=ytry,
  78.      psum:=psum+ptry-p[highest],
  79.      p:=replace(p,highest,ptry)) when ytry<y[highest]
  80.  
  81. -- shrink entire simplex
  82. shrinkp(k)[j] = .5*(p[k,j]+p[lowest,j]) dim[n]
  83. shrink = kk:=1,      
  84.   ((psum:=shrinkp(kk),
  85.     p:=replace(p,kk,psum)) when kk≠lowest,,
  86.     kk:=kk+1 ) while kk≤m,
  87.   y:=fp
  88.  
  89. -- define finite arrays so they can be assigned
  90. sump[j] = sum(p[ii,j],ii,1,m) dim[n]
  91. initp[i,j]= guess[j]+(deltas[j] when i=j,0) dim[m,n]
  92. fp[i]=f(p[i]) dim[m]
  93.  
  94. init = n:=count(guess), m:=n+1,
  95.        p:=initp,
  96.        psum:=sump,
  97.        y:=fp
  98.  
  99. replace(array,index,newval)[i] = newval when i=index,
  100.                                  array[i]
  101.  
  102. -- tuneable constants
  103. alpha=1; beta=.5; gamma=2
  104.